Breakdown of metastable step-flow growth on vicinal surfaces induced by nucleation 
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We consider the growth of a vicinal crystal surface in the presence of a step-edge barrier. For 
any value of the barrier strength, measured by the length £ ES , nucleation of islands on terraces is 
always able to destroy asymptotically step-flow growth. The breakdown of the metastable step- 
flow occurs through the formation of a mound of critical width proportional to L c ~ 1/V^es, the 
length associated to the linear instability of a high-symmetry surface. The time required for the 
destabilization grows exponentially with L c . Thermal detachment from steps or islands, or a steeper 
slope increase the instability time but do not modify the above picture, nor change L c significantly. 
Standard continuum theories cannot be used to evaluate the activation energy of the critical mound 
and the instability time. The dynamics of a mound can be described as a one dimensional random 
walk for its height k: attaining the critical height (i.e. the critical size) means that the probability 
to grow (k — » k + 1) becomes larger than the probability for the mound to shrink {k — > k — 1). 
Thermal detachment induces correlations in the random walk, otherwise absent. 

PACS numbers: 81.10.Aj, 68.55.Ac ,05.70.Ln 



I. INTRODUCTION 

The growth process of a crystal surface by deposition 
of particles resembles in many respects a phase separa- 
tion process. The role of the order parameter is played by 
the local slope u, whose average value m is fixed by the 
orientation of the substrate where deposition occurs. A 
uniform orientation is often not stable for a growing crys- 
tal, for thermodynamic^ as well as for kinetic? reasons: 
As deposition proceeds the flat profile gives way to a 
rough surface, with the formation and growth of mounds. 
The first stages of the destabilization process may occur 
through the appearance of large wavelength undulations 
or the formation of localized perturbations: the former 
happens if the orientation m is linearly unstable, the lat- 
ter if it is metastable. By varying m it is therefore pos- 
sible to pass from one regime to another, much in the 
same way as it is possible to go, in the phase separation 
of a binary alloy, from spinodal decomposition to droplet 
nucleation by changing the relative concentration. 3 

In the context of crystal growth, the unstable regime 
corresponds to a high-symmetry orientation, the counter- 
part of a symmetric quenching for a binary alloy. This 
regime has been studied in dozens of papers^ and sim- 
ilarities and differences with phase separation processes 
have been highlighted. In the following we study the so 
called vicinal regime^ when the slope m is sufficiently far 
from a high symmetry orientation: in this case the sur- 



face is made up of small terraces separated by steps and 
most of the deposited atoms attach to preexistent steps 
(step- flow growth) . This regime is the counterpart of the 
case where, for a binary mixture, the homogeneous phase 
is metastable, but it can be disrupted by the nucleation 
and growth of droplets of the stable phase. 

In the system considered in this paper, the key atom- 
istic ingredient of the destabilization is the presence of 
an energetic step-edge or Ehrlich-Schwoebel (ES) bar- 
rier— that hinders the descent of atoms from upper to 
lower terraces. This instability has a purely kinetic ori- 
gin, due to the nonequilibrium character of the growth 
process. The analytical investigation of a growing vici- 
nal profile reveals two possible linear instabilities: step- 
bunching 7 (the uniform density of steps is unstable) and 
step- meandering^ (straight steps are unstable). Step- 
edge barriers make the surface stable with respect to 
the bunching and unstable with respect to the mean- 
dering. This picture does not exhaust all the possibili- 
ties, because mound formation may occur also via non- 
linear mechanisms, leading to a late stage morphology 
which does not differ from what is obtained by growing 
on a (linearly unstable) high symmetry orientation. Two 
paths leading to mounding can be anticipated: atom- 
istic nucleation on terraces and defect formation, due 
to nonlinear meandering— In the following we consider 
the first possibility, which is the natural mechanism for 
two reasons: in the first place, it is solely responsible 
for destabilization if steps keep straight; 10 secondly, it is 
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the correct process to consider if we want to study the 
crossover from the unstable to the metastable regime, 
because atomistic nucleation plus step-edge barriers are 
responsible for destabilizing both a high symmetry and 
a vicinal surface. 

Since our interest is the destabilization via atomistic 
nucleation, we will consider a one- dimensional surface. 
This model may apply to a two-dimensional surface when 
steps are straight (see the final section for additional re- 
marks) . 

Let us now summarize the main features of the desta- 
bilization mechanism and our principal results. Freshly 
deposited atoms (adatoms) diffuse on terraces, generally 
until they are incorporated at steps. From time to time 
two of them meet and an atomistic nucleation occurs. 
Since the vicinal flat profile is linearly stable, the nu- 
cleated dimer is generally swept by the advancing steps. 
Only once in a while another nucleation event takes place 
on top of the dimer before it is reabsorbed, generating a 
very small mound. A mound may grow larger, but in 
general it tends to disappear, unless, by a fluctuation, it 
reaches a critical size such that it becomes more likely 
for it to grow rather than to shrink. At this point the 
flat profile is destabilized and it is not recovered by the 
dynamics. A crucial role is hence played by the critical 
nucleus, the localized perturbation of the profile (mound) 
that separates small perturbations that tend to be reab- 
sorbed from large ones that grow irreversibly. 

Section [H] highlights the inadequacy of atomistic ap- 
proaches based on the assumption that the critical nu- 
cleus has a fixed size, independent of the ES barrier. Sec- 
tion lllll shows that also the application of traditional con- 
tinuum methods for studying the breakdown of metasta- 
bility fails in our case. In Section llVl we discuss the form 
of the unstable current, which allows to distinguish be- 
tween the unstable and the metastable regime and gives 
insight into the crossover connecting them. 

In Sees. IV Al and IV Bl we present the results of Ki- 
netic Monte Carlo simulations showing that the surface 
is metastable for any strength of the additional step-edge 
barrier and that the instability time diverges exponen- 
tially for decreasing barrier. We also show that when 
the slope is sufficiently small we have a crossover to the 
linearly unstable regime, where the instability time di- 
verges as a power-law with the barrier. Sees. IV CI and 
IV Dl present an analysis of the properties of the critical 
nucleus, whose width is related to the instability of a 
singular surface, along with a detailed study of how the 
nucleus is dynamically generated, allowing a quantitative 
derivation of the instability time. 

Finally, in Sec. IVII the global physical picture is 
shown to be robust when thermal detachment is allowed: 
metastability holds for any barrier strength also in this 
case. Nonetheless, we show that thermal detachment has 
nontrivial effects on the development of the critical nu- 
cleus and this results in a remarkable increase of the 
instability time. A discussion (Sec. IVIIfl concludes the 
paper. 



In a recent Letter^ii we have presented some partial 
results of this line of research. Though some overlap is 
unavoidable, we have minimized it in the present paper 
by referring to the Letter when needed and focusing here 
on completely new results: the failure of the continuum 
approach fSec. Illl(l : the study of the crossover between 
singular and vicinal regimes fSec. ITVjl : the detailed nu- 
merical and analytical investigation of the destabiliza- 
tion process at the atomistic level (Sees. IV CIV Dl and 
the Appendix); the simulations with thermal detachment 
(Sec. ED . 



II. THEORY WITH A CRITICAL NUCLEUS OF 
FIXED SIZE 

A simple and appealing argument for analyzing the 
destabilization of step-flow on vicinal surfaces has been 
recently introduced by Kallunki and Krug2*A£ and it is 
based on the idea that the critical nucleus is a mound 
of height two. The stability or instability of step-flow is 
therefore assessed by evaluating whether such a nucleus 
is formed or not by the dynamics. This is performed 
by comparing the mean timescales of the two processes 
that can occur once an island is formed on a vicinal ter- 
race: either the island is reabsorbed by the advancing 
step (and this requires a time i adv ) or a second-layer nu- 
cleation event takes place on top of it (this occurs over a 
time t 2 nd)- If ^2nd < ^adv a mound of height two is formed 
and this leads to the destabilization of the surface. 

There are several ways of implementing in practice this 
criterion and they all lead to results that are equivalent, 
apart from small differences in the numerical factors. We 
refer here to the simplest model of epitaxial growth, with 
deposition occurring at rate F, adatom diffusion at rate 
D and the strength of the ES barrier quantified by the 
length ^ ES = D/D' — l, where D' is the interlayer diffusion 
coefficient. One way to implement the criterion is to 
impose that the mean number of second-layer nucleation 
events occurring during the deposition of a monolayer is 
larger than 1. The number of such nucleation events can 
be evaluated by integrating the second-layer nucleation 
rate uj over time up to l/F, 

f l/F 

N nuc (t = l/F) = / dsu[R(s)], (1) 
Jo 

where R(s) is the size of the island at time s and 
p2n4 / p \ 

^> = «w(i+6f)> < 2 » 

(k being a constant of the order of oneji&ii that will be 
neglected in the following). 

Assuming^ that half of the deposited adatoms con- 
tribute to the growth of the island (and the other half 
is incorporated at the advancing steps), we have R(s) = 
F£s/2, where I = 1/m is the size of the vicinal terrace. 
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Performing the integrals and imposing N nuc (t — l/F) > 
1, we obtain 

4s > l c ES = M^f - ±£. (3) 

We stress again that other implementations of the same 
criterion yield virtually identical results. For instance 
one can impose that the probability p nuc (t*) that second- 
layer nucleation occurs within a time t* = l/(2F) after 
the creation of the island is larger than 1/2. In this way 
we obtain the same condition J2J except for a factor ln(2) 
multiplying the first term on the r.h.s.. 

We conclude that the assumption of a critical nucleus 
of height two implies the existence of a threshold: for 
small enough barrier, step-flow growth should be fully 
stable with respect to the nucleation of mounds. As 
shown below, this prediction is in striking contrast with 
numerical simulations, showing that, even for values of 
£ ES three orders of magnitude smaller than £ ES , the sur- 
face is in fact metastable. This is a clear indication that 
something in the argument just presented is not correct. 

One could think that a higher critical nucleus would fix 
the problem. However, a critical nucleus of fixed height 
k c would always imply the existence of a range of small 
values of £ ES for which step-flow growth is fully stable. 
This can be seen by considering that, for any k c , the con- 
dition for instability is of the form A(£ ES ) > Ao, where 
A is a growing function of £ ES , expressing how likely the 
formation of the critical nucleus is, and Ao is a constant, 
whose precise value depends somewhat arbitrarily on the 
details of the implementation. If A(£ ES — 0) > Ao then 
the growth would be unstable even for £ ES — and this 
is absurd. Hence A(£ E3 — 0) must be smaller than Ao 
and this implies that in an interval of small £ ES growth 
is stable. Metastability down to £ ES — would be pos- 
sible only if A(£ ES = 0) exactly equals Ao, but this is 
impossible, given that a precise value of Ao cannot be 
unambiguously defined. 

Another reason for the failure of the atomistic ar- 
gument discussed above could be that comparing the 
average values of t 2ad and i adv is not appropriate, we 
should consider instead their full distributions Pi(t adv ) 
and P2(i 2 nd)- The surface is unstable if there is a finite 
probability that t 2nd < i adv , i.e. if 

Pt2»d<tadv = / dt 2nd P 2 (t 2nd ) di adv Pl(t adv ) >0. 

Jo Jt 2nd 

(4) 

Let us now consider the limit £ ES — ► 0: both distribu- 
tions have a well defined limit, because the nucleation 
rate on top of a terrace and the velocity of a step reach 
finite values for £ ES — 0. Since P2(t 2nd ) > for small 
but finite t 2ad , the quantity Pt 2nd <t adv is always greater 
than zero. In simple terms, considering the condition (0J 
for metastability leads to the wrong conclusion that the 
vicinal surface should be destabilized for £ ES — as well. 
In conclusion, the hypothesis of a critical nucleus with 
fixed height is not compatible with numerical results and 



known properties for £ ES — 0. As it will be shown below, 
the diverging size of the critical nucleus is the crucial el- 
ement of the breakdown of metastable step-flow growth. 



III. CONTINUUM CAHN-HILLIARD THEORY 
FOR THE METASTABILITY 

The goal of this section is to study the metastability 
by applying the Cahn-Hilliard (CH) theory developed for 
ordinary phase-separation to a minimal continuum model 
for epitaxial growth. We will show that this approach is 
completely unable to explain the dynamical behavior of 
the growing surface even in qualitative terms. A very 
short introduction to the continuum approach is given in 
the next two paragraphs. 

In the absence of events which do not conserve matter 
(evaporation) or volume (formation of overhangs), the 
dynamics of a one-dimensional growing crystal surface 
is described by dtz = —d x J, where J is the current de- 
scribing all microscopic processes occurring at the surface 
and z is the local height h with respect to a comoving 
frame, z = h — Ft (F is the flux intensity). Since we 
are not studying phenomena depending on the substrate- 
adsorbate interaction, J can be safely assumed to depend 

only on spatial derivatives of z, U = z', u' , u", In the 

next section we will evaluate numerically the slope depen- 
dent current, j ES (u), also called step-edge or ES current 
because it is due to asymmetric sticking to steps. j ES (u) 
is linear at small sloped attains a maximum for u = mo 
and finally decreases. In addition, other terms depending 
on higher order derivatives contribute to the current J. 
The most important has the form Ku", is called Mullins 
term, and is due to several microscopic processes^ ther- 
mal detachment from steps and the random character of 
nucleation and diffusion. In the following we will consider 
a minimal model, with J = Ku" + aj ES (u) . The prefac- 
tor a is chosen so that j' ES (0) = 1: it might be included in 
the definition of j ES (u), but its explicit appearance will 
be useful later. 

Given the form of J, the evolution equation for the 
surface slope u has the form of a generalized CH equation 

d t u=-dl[Ku"(x)+aj ES {u)]. (5) 

Linear stability analysis of a vicinal surface of slope 
to, z = mx + 5zexp(u)t + iqx) provides the spectrum 
w (<?) = a 3 ES ( m )<l 2 - K( f- H m > TOq, j' ES {m) < and 
the surface is linearly stable; if m < mo the surface is lin- 
early unstable, the most unstable wavevector (the one for 
which uj(q) is maximum) is q c = y // aj ES (m)/2K , the most 
unstable length is L c = 2n/q c = 2n 2K / aj' ES (m) , and 
the linear instability emerges after a time t c i=s l/u)(q c ) 
K / (aj' ES (m)) 2 . For small barriers 2 a is proportional 
to £ ES ; therefore the instability appears with a typical 
lengthscale of order l/v / 4s an d after a typical time of 
order l/£ ES - These temporal and spatial scales also di- 
verge when m — > TOq , because j' ES (m) — > 0. 
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Stationary solutions of Eq. JSJ) are given by Ku" (x) + 
°ijEs( u ) — M> where \i is related to m by the condition 
(u) = m. The critical nucleus is a localized steady state, 
characterized by u(x) — m everywhere except in a region 
of finite size. This implies that, for the critical nucleus, 
^ = aj ES (m) and 

Ku"(x) + a[j ES (u) - j ES {m)} = 0. (6) 

Rather than reproducing the standard procedure for 
getting the energy AT of the critical nucleus^, we simply 
remark that rescaling x allows to find how AT depends 
on a and K. If X = \f%x, Eq. © has the parameter- 
free form uxx + \j-Es(u) — j ES (m)] = and the pseudo 
free-energy T [such that d t u = d\ (f^-)] is 

T = J dx[^{d x uf + aU(u)} 

= Vcdt J dX[^(d x u) 2 + U(u)} (7) 

where U'(u) = —j ES (u). In conclusion, 

AT{a,K) = VaKAT(l,l). (8) 

In a thermodynamic system, once the activation en- 
ergy AT has been calculated, the instability time is given 
by r ina ~ exp(AJ 7 /fcsT), because the system overcomes 
the barrier by thermal fluctuations. In the far-from- 
equilibrium case of a growing surface, fluctuations allow- 
ing to system to leave the local minimum have different 
origins: fluctuations in the flux (shot-noise) and fluctu- 
ations in the nucleation and diffusion processes. Their 
amplitude plays the role of ksT in the expression for 

Tina ■ 

The analytical form of j ES (u) , valid for any u and £ ES , 
is not easy to write, but we need here only the limiting 
form for vanishing barriers and finite slope. With these 
caveats, all the £ ES dependence is contained in the pref- 
actor a, which is proportional to £ ES . Hence for fixed 
slope, the unstable current is linear in the ES length in 
the limit of small barrier. According to Eq. (JSJ), the acti- 
vation barrier is therefore an increasing function of ^ ES , 
the opposite of what is expected and observed! This sur- 
prising result cannot be changed by taking into account 
the possible effects of £ ES on the noise amplitude, because 
nucleation noise has a finite limit for £ ES — 0. We must 
conclude that the continuum minimal model, Eq. (JSJ), is 
not appropriate for the study of the breakdown of step- 
flow growth. Additional comments are deferred to the 
Conclusions. 



IV. CROSSOVER FROM SINGULAR TO 
VICINAL FOR A TILTED SURFACE 

Before reporting the results of a numerical study of the 
breakdown of step-flow growth, let us discuss what occurs 



at the boundary between singular and vicinal surface ori- 
entations. Linear stability analysis predicts two distinct 
types of dynamical behavior, depending on how large is 
the average global tilt m with respect to the value too for 
which the Ehrlich-Schwoebel nonequilibrium current j ES 
is maximal. For to < too (singular surfaces) a flat profile 
is linearly unstable for q < q c (see the previous section). 
For to > Too instead (vicinal surfaces), perturbations of 
small amplitude decay, since the flat profile is linearly 
stable, and, if any destabilization occurs, this has to be 
triggered by fluctuations of sufficiently large amplitude. 

In order to study the way the destabilization of vici- 
nal surfaces takes place, it is therefore crucial to check 
that the surface considered is effectively vicinal, i.e., its 
slope is larger than too. This point is made complicated 
by the fact that Too may in principle (and it actually 
does) depend on £ ES , so that it may happen that a sur- 
face with fixed m is singular or vicinal depending on the 
value of £ ES . The variation of too is not expected to be 
large, because the form of the current predicts that the 
ES current is maximal for finite slopes of the order of 
1 /£ D (see Ref. @) for both £ ES — * and £ ES — ► oo (toq 
and TO^, respectively). In any case, this effect must be 
taken into account for a detailed understanding of the 
destabilization process. For this reason we have numer- 
ically evaluated the form of the nonequilibrium current 
j ES (u), by measuring, for the KMC model described be- 
low, the imbalance in the number of hops toward left or 
right that adatoms take during the growth of the first 
few monolayers. 3,7 More precisely, if Ni and N r are the 
number of steps toward left and right, respectively, then 
the current is 

\Ni-N r \ 

J ES = 77 , (9) 

where N a d is the number of adatoms deposited. 

The results, reported in Fig. ^ display a very slow 
crossover between too ~ 0.07 for small £ ES and to « 0.03 
for the largest value of £ ES we could consider: we stress 
that too varies by a factor two against a variation of £ ES 
over five orders of magnitude. As expected, these values 
are of the order of 1/£ D w 0.025. 

We can conclude that three types of behavior exist for 
the system under consideration (Fig [2J) • For to < to x> « 
0.03 the surface is singular (hence linearly unstable) for 
all values of the ES barrier. For to > to[J ~ 0.07 the 
surface is always linearly stable and metastability can be 
studied down to ^ ES — > 0. For intermediate slopes the 
surface is vicinal for large £ ES , while it becomes linearly 
unstable as ^ ES is decreased. In this case a crossover is 
expected in the way the instability time depends on £ ES . 

Notice that the value m = 1/15 ~ 0.067 considered in 
the following sections just falls within this interval. It is 
important to remark (Sec. IHIfl that the linear instability 
for a singular surface does not imply that the instability 
develops immediately. The time required for a sizeable 
amplification of the initial fluctuations depends on to and 
it diverges for m — > toq or £ ES — > 0. 
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FIG. 1: Main: Plot of the nonequilibrium current j'es(w) 
as a function of the slope m, for several values of the 
Ehrlich-Schwoebel length. From bottom to top: l ES = 
0.1, 0.4, 4, 40, 400, 4000, 40000. Curves have been shifted along 
the vertical axis for clarity. Inset: Plot of the value mo of the 
slope for which the current is maximal. 
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FIG. 2: Schematic plot of the different regimes as a function 
of the slope m and of the ES barrier. 



FIG. 3: Double logarithmic plot of ln(r ina ) as a function of 
^es, for several values of the slope m and of the thermal de- 
tachment rate Rd- Error bars are smaller than symbol size. 



neighbor belongs to the same terrace of the first site, at 
rate D 1 < D (because of the Ehrlich-Schwoebel effect) if 
the atom must descend a step. Dimers and larger islands 
are immobile. Thermal detachment is possible at rate 
Rd- When a particle is detached from a step or an island 
it is placed on the lower terrace at distance 1 from the 
step. We have checked that if detachment events into 
the upper terrace are also allowed (with the same rate 
of the detachments into the lower terrace), results do not 
change appreciably. The dynamics has been implemented 
via a rejection-free type of algorithmpiSiiS 

We have always kept fixed the deposition rate F = 1 
and the diffusion coefficient D = 5 x 10 5 , so that the 
diffusion length is £ D w 40, as in Ref. [2(]. The size of the 
system has been always L — 1000, large enough (L 3> A c , 
see below) to allow for unstable mounds to form. We 
have varied D', in order to change the ES length ^ BS = 
D/D' - 1, and R d . 



B. The instability time 



V. KINETIC MONTE CARLO SIMULATIONS 

A. The model 

We perform Kinetic Monte Carlo (KMC) simulations 
of the simplest model for epitaxial growth in d = 1 . The 
surface is represented by a set of integer height variables 
hi, defined on a one-dimensional lattice of L sites. An 
average tilt to is imposed through helical boundary con- 
ditions, hi+L — hi + mL. The average terrace size is 
thus £ = 1/m. The initial condition is a regular train 
of steps, hi — [|] , where [x] is the integer part of x. 
Deposition (hi — ► hi + 1) occurs at rate F on randomly 
selected sites. Singly bonded atoms (adatoms) attempt 
diffusion hops to nearest neighbor sites, at rate D if the 



In order to analyze the breakdown of metastability, 
we focus first on the destabilization time r ins , i.e. the 
time needed for the formation of the first critical nucleus. 
The identification of this object is not trivial and it has 
been discussed in Ref. [llj. Using the same numerical 
procedure for the identification of the nucleus, we have 
computed r ins for several values of the slope to, of the 
rate for thermal detachment Rd, and with varying t. E g . 

The double-logarithmic plot of ln(r ins ) vs £ ES (Fig. 01 
indicates a divergence as £ ES — > 0, with no clear-cut expo- 
nential or power-law behavior. Common (and intuitively 
expected) features are that increasing to induces a de- 
lay in the destabilization process, just as, for fixed to, 
thermal detachment does. 

More insight into how the destabilization process de- 
pends on £ ES is obtained by considering the formation 



6 




FIG. 4: Main: Double logarithmic plot of ln(l/p„) as a func- 
tion of £ E s , for several values of the slope m and of the thermal 
detachment rate Rd- Error bars are smaller than symbol size. 
A straight line indicates a divergence as exp(a/vSs)- In- 
set: 1/pu for m = 0.067 and Rd = is plotted in a double 
logarithmic-plot, showing that the divergence for £ ES — > is 
close to the l/£ls behavior (solid line) predicted for a singular 
surface. 

of the critical nucleus as the combination of two distinct 
processes. In the first place a dimer must be nucleated 
on a vicinal terrace. Then the dimer must grow up to a 
mound of a certain critical size that makes it unstable. 
The first event occurs at a rate l/r nuc i(L), which is pro- 
portional to the system size L and is independent of the 
ES length in the limit £ ES — ► 0. The second event occurs 
only with a small probability p u , strongly dependent on 
£ ES but not on L. Hence we can write 

r ma (i,4s) = T ona (L,£ ES )—^—. (10) 

Pu{£ ES ) 

The value of r nucl can be written in its turn as 
Tdim(^, £ E s)£/L where r dim is the time to nucleate a dimer 
on a vicinal terrace of size £ and L/£ is the number of such 
terraces: T dim could be computed analytically following 
Ref.Q but it is straightforward to evaluate it numeri- 
cally. Furthermore, an accurate determination of t„ uc1 is 
important in order to evaluate correctly 1 / p u (£ ES ) , which 
can be extracted by plotting r ins /T nucl vs £ ES (Fig. 0}. 

It turns out that, in almost all cases, l/p u diverges 
with £ ES — + as exp(a/£2 s ), with 7 w 0.5 and without 
evidence of a finite threshold. The only exception is the 
case of weaker slope, to = 0.067, where the divergence for 
small barrier is slower than exponential. This different 
behavior can be understood on the basis of the form of 
the unstable current j ES reported in Fig. ^ For ^ ES < 10 
the maximum too of the current is very close to the slope 
to = 0.067 considered here. As a consequence, as the 
barrier is reduced the surface is more and more singular, 
rather than vicinal. The destabilization takes place be- 
cause of a (weak) linear instability instead than through 
the breakdown of a metastable state. As a consequence 
the dependence of 1/p u on £ E g is not exponential, cis for 
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FIG. 5: Plot of the width A c of the critical mound as a func- 
tion of the ES length for several values of the slope m and 
the thermal detachment rate Rd- Error bars are smaller than 
symbol size. The full circles are the values of the length L c 
associated to the linear instability for a singular surface. 

the other curves in Fig. 0] The divergence is a power- 
law, with an effective exponent 2.2 (see Fig. QJ inset), 
not far from the behavior 1 / £^ s expected from the linear 
instability of a singular surface. 

C. The mound width and its critical value 

The evolution that, starting from a dimer nucleated 
on a vicinal terrace, leads to the formation of an unsta- 
ble mound, is a very complicated process depending at 
each time on the detailed form of the profile. In order 
to understand it theoretically, one must focus on few pa- 
rameters. The features that most naturally characterize 
a mound are its height k and its width A. Both are use- 
fid in the analysis of the destabilization process. We first 
discuss the role of the mound width, which, despite be- 
ing numerically harder to determine, it is conceptually 
more fundamental. In the following subsection we will 
consider the temporal evolution of the height, which is 
easily measured in simulations. 

As discussed in Ref. 0, a good definition of the mound 
width is the size of the terrace immediately below the top 
one. For each value of to, £ es and Rd it is thus possible 
to measure the critical width A c , as the average width of 
mounds that start to diverge irreversibly. A plot of this 
quantity (Fig. [5J reveals the remarkable property that 
A c does not depend on to, nor on the rate of thermal 
detachment Rd- Moreover the dependence on the ES 
length is a power-law, with an effective exponent not far 
from 1/2. We now focus on the case without thermal 
detachment leaving for the next section the analysis of 
the effect of a finite Rd- 

For the case without thermal detachment, the behav- 
ior displayed in Fig. [S] is naturally and consistently in- 
terpreted^ by relating the critical width A C (£ ES ) to the 
critical size L C (£ ES ) associated to the linear instability of 
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a singular surface 



A c (4s) oc L c {i 

E 



(11) 



This interpretation is further corroborated by a direct 
numerical evaluation of L C (£ ES ) as the smallest system 
size such that for a singular surface the instability can 
take place. The comparison of A c with L C (£ ES ) (Fig. 
strongly supports the validity of Eq. , confirming that 
a growing mound becomes unstable when its width is 
large enough to trigger the same type of instability at 
work for a singular surface. In this respect the width A is 
the fundamental quantity determining whether a mound 
is stable or unstable. Because of the difficulties in its 
numerical determinationii, it is not easy to analyze the 
detailed evolution of A in order to understand how the 
critical nucleus is generated by the dynamics. In partic- 
ular, we aim at deriving a quantitative formula for the 
instability time T ina or the relevant factor l/p u - For this 
purpose it is more useful to study the mound height fc, 
which is easily measured in simulations. 



D. Evolution of the mound height 

The height is not as fundamental as the width for 
the determination of whether a mound is supercritical 
or not. For example a critical height cannot in princi- 
ple be rigorously defined. This can be seen by taking a 
mound generated by the dynamics just a bit wider than 
A c and strongly reducing its height: it remains super- 
critical. Nevertheless the dynamics couples the mound 
width and height so that the evolution of k tipically re- 
flects that of A and a height k c of the critical nucleus can 
be in practice defined. 

If a mound is described only in terms of its height, its 
dynamics is reduced to a one-dimensional random walk 
(RW), with only two possible processes for a mound of 
height k: it may either become taller (k — * k + 1) via 
dimer nucleation on its top terrace, or shorter (k — > k— 1) 
if the bottom terrace delimiting the mound is filled. The 
increase (decrease) happens with probability p+(k) [1 — 
p+(k)}. It is possible to determine, via KMC simulations, 
the function p+{k), which is reported in Fig. [HJ 

The stochastic nature of destabilization for large k is 
evident, consequence of the shape of p+(k). For small k, 
the probability that the mound grows is much smaller 
than 1/2. This implies that only rare events lead to 
mound growth. This effect is reduced as k is increased, up 
to a certain value k c , such that p+(k) > 1/2 for k > k c . 
The threshold k c plays the role of the effective height of 
the critical nucleus. The evolution is therefore biased to- 
ward the flat profile for k < k c while the bias is toward 
higher k (destabilization) for k > k c . This translates 
into microscopic terms the metastable nature of a vici- 
nal surface: it is stable with respect to small fluctuations, 
but unstable when a large fluctuation by chance appears. 
Fig. shows that k c grows as £ ES is reduced, invalidating 
the assumption of the argument in Sec. ITU 
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FIG. 6: The probability p+(k) that a mound of height k even- 
tually turns into a mound of height k+ 1, rather than decay- 
ing to a mound of height k — 1. Results are for m = 0.2 and 
Rd = 0- 



Using the values of p+{k) extracted from the full KMC 
simulations, and taking them as the transition probabil- 
ities for an uncorrelated one-dimensional random walk, 
it is possible to compute l/p u as the fraction of times a 
walker starting in k = 1 reaches a large reference value 
(larger than k c , say 15) without touching the absorb- 
ing boundary k — 0. The results, presented in Fig. 
agree very well with the outcome of simulations, indi- 
cating that, despite its rather crude simplifications, the 
random walk picture captures much of the destabiliza- 
tion process, and it allows to recover the dependence of 
the instability time r ins on £ ES . It is important to re- 
mark that this result is not at all trivial. Correlations 
in the random walk can generate, with the same tran- 
sition probabilities p+(k), very different values of l/p u . 
The agreement exhibited in Fig. means that the un- 
correlation of RW steps is a correct assumption. As it 
will be shown below, this is no more true when thermal 
detachment is allowed. 

Within the RW picture it is also possible to derive an- 
alytically a formula for l/p u . As shown in the Appendix 
one has to solve a stationary convection-diffusion equa- 
tion with a space-dependent drift related to the form 
of p+(k). By looking at Fig. it is reasonable to take 
the explicit espression p+{k) = poe k for k < k c and 
p+{k) — p max for k > k c , with poe kc = p max . In this 
way we obtain (see the Appendix), in the limit k c 3> 1 



1 



2p max -l 

(e 2 -l) 



(12) 



The behavior of 1 /p u with decreasing £ ES is then due to a 
power-law divergence of the critical height k c . Since the 
height and the width of the critical nucleus scale in the 
same way2i it is then possible to write l/p u ~ e a ( m ) A <=. 



8 




FIG. 8: The probability p+(k) that a mound of height k even- 
FIG. 7: Comparison between the values of l/p u obtained in tuflUy tumg into a mound of height k + ^ rather than decay _ 

numerical simulations (empty symbols) and the same quantity ing to a mound of height k _ 1 Results are for m = . 2 and 

computed within the uncorrelated random walk with transi- _ cjqq 

tion probabilities p+(k) (filled symbols). 



VI. THE EFFECT OF THERMAL 
DETACHMENT 

We now analyize how the previous picture of the break- 
down of the metastable state is modified by the possibil- 
ity that atoms detach from steps or islands with rate Rd- 
Figures |3| EI and [5] show the effect of Rd on the insta- 
bility time and on the critical width A c , respectively. At 
first sight these results might seem contradictory. While 
the critical nucleus width does not change appreciably 
when Rd grows from to 70 (Fig.JSJ, l/p u increases by 
a factor larger than 10 3 (Fig. ^J. 

In Sec. IV CI we related A c to L c , the minimal length 
for the onset of the (linear) instability on a singular sur- 
face. A numerical check shows that also L C (£ ES ) does 
not change when thermal detachment is allowed: we can 
therefore conclude that Eq. remains valid also in the 
presence of thermal detachment. From the continuum 
theory (Sec. Illlf) . L c ~ ^K/l ES , where the coefficient 
K contains several contributions, one of which is due to 
thermal detachment, K = Kq + K t her- The fact that 
L c is independent of Rd implies that for Rd < 70 the 
thermal contribution K t her is negligible. 

To clarify why thermal detachment has instead dra- 
matic effects on l/p u , we have computed, exactly as for 
Rd = 0, the transition probabilities p+{k) for the sim- 
plified random walk picture (Fig. |3J). The presence of 
thermal detachment leaves the shape of p+{k) qualita- 
tively the same. The value of the threshold k c (the crit- 
ical nucleus height) does not change much with Rd, in 
agreement with the result A c : the critical nucleus does 
not practically depend on Rd- However, for the same k, 
p+(k) is reduced as Rd grows, in particular for small k. 
This tendency has a simple interpretation in terms of a 
microscopic process induced by thermal detachment: the 
decay of the freshly nucleated dimers on top of mounds. 
While for Rd = a dimer on top of a mound is stable 



and can only grow by incorporating incoming adatoms, 
for Rd > a dimer can split in two adatoms, which 
may diffuse and be incorporated by surrounding steps, 
restoring the situation with no dimer. This decay re- 
duces the probability p+{k) that a mound grows higher. 
A smaller p+(k) clearly implies an increased l/p u . Con- 
versely, when a dimer splits during the growth of a high- 
symmetry surface, adatoms are not lost because there are 
no preexistent steps which capture them: hence such a 
splitting has no dramatic effect, leaving L c (and hence 
on A c ) unchanged. This is the reason why switching Rd 
on is relevant for T ina and not for A c . 

The reduction of p+(k) shown in Fig. [H]does not tell 
the whole story about the effect of thermal detachment. 
If we compare the value of l/p u determined numerically 
for the uncorrelated random walk with transition prob- 
abilities p+(k) with the same quantity measured from 
the full KMC simulations (Fig. 0), we realize that the 
results do not agree and the mismatch grows with Rd- 
The smaller p+{k) is not enough to explain the instabil- 
ity time measured in KMC simulations: the assumption 
of uncorrelated steps is wrong when Rd ^ 0. The origin 
of the correlation is another subtle but relevant effect of 
the decay of dimers. When a dimer is formed on the top 
terrace of a mound (k — > k + 1), owing to thermal de- 
tachment it is very easy that the dimer will soon decay 
so that k + 1 goes back k. This increased chance that a 
growth event is immediately followed by a decrease of k 
invalidates the assumption of uncorrelated steps. There- 
fore, the uncorrelated random walk picture is no more 
useful for determining quantitatively l/p u and from it 
the instability time. 



VII. CONCLUSIONS 

In a nutshell, we have investigated the breakdown of 
the metastable step-flow for epitaxially grown vicinal sur- 
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faces. We have performed Kinetic Monte Carlo simula- 
tions showing that metastability holds for any strength of 
the ES barrier: the time needed for destabilization to oc- 
cur diverges exponentially when the barrier goes to zero, 
both in the presence and in the absence of thermal de- 
tachment. Previous atomistic and continuum approaches 
to the problem do not reproduce this phenomenology. 
The crucial event determining the end of the metastable 
state is the formation, via a rare fluctuation, of a critical 
nucleus whose size diverges as the barrier goes to zero. 
The dependence of the instability time on the parameters 
of the system can be summarized as 



,(L,m,£ ES ,R d ) 



,{m,£ ES ,R d ) [a(m,R d )L c (e ES )] 



(mL) 



. ( . 13) 

where r dim is the time to nucleate a dimer on a vicinal 
terrace and L C (£ ES ) is the length associated to the linear 
instability on a singular surface. r dim has a well-known 

analytical form— (at least for R d = 0) and L c = yj 

The only quantity which has not been determined is the 
prefactor a in the exponent. The £ ES -dependence in T dim 
is weak and negligible for small £ ES . 

The results presented fully clarify how the destabiliza- 
tion of step-flow comes about and also provide a semi- 
quantitative picture. The analytical interpretation would 
be completed by a direct computation of the probabilities 
p+(k) based on the microscopic dynamics. This remains 
a challenge for future work. 

Another intriguing open question has to do with the 
difficulties of the continuum approach. In Sec. lIIII we have 
shown that the minimal model with J — Ku xx + aj ES (u) 
is completely unable to explain the phenomenology of the 
metastability process. A likely explanation for the failure 
of the continuum theory is that the current J considered 
is too simple: it is well known that terms breaking the 
z — > — z symmetry may be relevant, and a look at pro- 
files from KMC simulations indicate that they play some 
role. These symmetry-breaking terms usually make the 
equation of motion nonderivable from a potential so that 
energetic arguments as those presented in Sec. IIIII are 
not applicable. An analytical approach to the metasta- 
bility of generalized Cahn-Hilliard equations with sym- 
metry breaking terms is still to be found. But solving 
this problem may not be enough. Because of such addi- 
tional terms, non-analytical regions could appear in the 
surface profile, in particular in the surroundings of min- 
ima where mounds merge. If this occurs a continuum 
approach may not be altogether applicable. 

Despite these possible objections to the use of the mini- 
mal model © it is important to remind that such a model 
is instead appropriate to describe at least qualitatively 
the linear and the nonlinear regimes of the destabiliza- 
tion of a singular surface. Why the model is not appro- 
priate for the understanding the decay of the metastable 
regime is, in this light, a puzzle. 

Finally, some considerations on two-dimensional sys- 
tems. The argument given in Sec. [HJ based on the hy- 



pothesis that critical mounds have constant height, can 
be easily extended to d — 2 and still provides a thresh- 
old. Within our approach, we can write the analogue of 
Eq. JTDJ, 



r iM (L||,Lx,4s) 



L ± (L«/£)p u ' 



(14) 



where Lii and L± are the linear size of the system in the 
direction parallel and perpendicular to the slope, respec- 
tively. T dim (£) is the time to form a dimer on a terrace 
of size £ in the parallel direction per unit length in the 
perpendicular direction. Again, p u {£es) is the probabil- 
ity that a dimer grows up to the critical size instead of 
being absorbed. Assuming that the form of the probabil- 
ities p+(k) will be similar to the one-dimensional case, it 
is reasonable to expect l/p u ~ exp(fc c ) also in this case. 
Since k c scales as A c in d = 2 as well, we conjecture that 
also in two dimensions the instability time diverges as 

1/2 

exp(a/£ ES ) when the ES barrier goes to zero. 



APPENDIX A: RANDOM WALK MODEL 

In this Appendix we compute analytically the prob- 
ability p u that a freshly created dimer evolves into an 
unstable mound instead of being reabsorbed, given the 
probabilities p+(k) [1 — p + (k)] of the process k — > k + 1 
[k — > k — 1] . The problem is equivalent to the compu- 
tation, for a random walk with absorbing boundaries in 
and N — ► oo, of the exit probability in N, with the 
walker starting initially in k = 1. 

Let us call e(x) the exit probability if the walker is 
deposited in x. This quantity obeys the recurrence rela- 



tio: 



,22 



e(x)=p + (x)e(x+l) + (l-p+(x))e(x-l), (Al) 

with boundary conditions e(0) = and e(N — > oo) = 1. 
In the continuum limit 



= a(x)e'(x) + l -e"{x), 



(A2) 



with a(x) — 2p + (x) — 1. Therefore we must solve 
a stationary convection-diffusion equation with space- 
dependent drift. The equation can be solved by sepa- 
ration of variables yielding 



i(x) = Ci +C / dx'e- 2f ^'K 



(A3) 



with f(x') = f* a(x")dx". 

To proceed further we specify the form of p + (k), in 
agreement with Fig. 

= { P p1 tit (A4) 
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with poe kc — p max . Using this form we have, for x < k c , 

f(x) = fx(x) = 2p e x - x, (A5) 

so that, for x < k c , 

e(x) = e 1 (x) (A6) 
(1 + 4p )e 4po - (1 + 4poe x )e- 4 P° e!C 



= Ci + Co 
while for x > k c 

e(x) = e 2 (x) = ei(k c ) + C Q - 



16pg 



-2/i(fe c ) 



(A7) 



2(2p max -l)" 

The boundary condition e(0) = sets C\ = 0. The 



second boundary condition determines the value of Co 



_L_ _ (1 + 4p )e 4 P" - (1 + 4 Po e x )e- 4poe 
Co ~ " 



16^ 



e -2/i(fec) 

"2(2p m „-l)' 
(A8) 

The quantity we are interested in is l/p u = l/ei{x — 

1). In the limit k c ^> 1, which implies po — * 0, and 

e -2/i(fe c ) ^ e 2fc c ^ we g n j 



1 



1 + 



->/<„ 



(e 2 -l) 



„2fe c 



(A9) 
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